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PARAMETRIC STUDY OF THE DYNAMIC STABILITY OF TOWED VESSELS 



Several accidents in towing operations of barges or 
disabled ships in restricted and open waters have made 
necessary the investigation of the course keeping stability 
of towed vessels. In this work a non-linear mathematical 
model is used to simulate the slow surge, sway, and yaw 
motions of a vessel towed by a heavy catenary towline. The 
effect of geometric parameters of the system on the 
stability of equilibrium configurations is analyzed. 

It is shown that for certain choices of towing system 
parameters, dynamic loss of stability may occur which 
results in qualitatively different asymptotic response. The 
results of this study identify regions in the parameter 
space that lead to either safe operations or hazardous 
system response. 
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I. INTRODUCTION 



A. BACKGROUND 

A long history of towing accidents resulting in loss of 
life, damage to property, and pollution of the environment 
have prompted many studies into the dynamics of towing 
operations. Of primary concern was the motions of the towed 
vessel in the horizontal plane (yaw, sway, and surge) . 
Excessive and unstable motions could lead to collisions and 
capsizing. The ability to predict the motion of a 
particular towing system would be of particular benefit to 
ship designers and towing operators, by identifying those 
situations where the towing operations would be the safest, 
or those which must be avoided. 

Previous studies at the University of Michigan and 
elsewhere had developed mathematical models and numerical 
techniques to analyze towing dynamics, and had identified 
those parameters which are of primary importance to the 
stability of the towing system. The linear model usually 
used to describe ship motions [Ref. 1, Chapter 7] is 
inadequate for the towing problem. Non-linear models, as in 
[Ref. 2], must be utilized to accurately describe the 
towing system. These studies had identified the position of 
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the towline attachment point on the towed vessel and the 
towline tension as the most significant (and controllable) 
parameters of the towing system. 

B . PROBLEM CONDITIONS 

In this study, computer programs developed in [Ref. 3] 
were used to analyze the effect of different parameter 
combinations on the towed stability of three vessels. These 
programs use bifurcation to identify the unstable and 
stable regions of the parameter space. The principal 
parameters studied were (Fig.l): 

1. longitudinal position of the towline attachment point 
forward of the towed vessel's center of gravity, xp; 

2. athwartships position of the towline attachment point 
port or starboard of the towed vessel's centerline, 

yp; 

3. length of the towline, Lw. 

In the model used in this study, unlike [Ref. 2], the 
towline is modeled as an inextensible caternary, thus 
making towline tension a function of its length. The model 
conditions were: 

1. speed of towing vessel of 2 knots; 

2. towing vessel on steady course; 

3. calm seas, no wind; i.e., no external environmental 
forces . 

Characteristics of the towed vessel were inputted into the 
programs from a data file containing hydrodynamic 
coefficients, resistance data, towline characteristics, and 
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> Figure 1. Problem Geometry 
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skeg data, if applicable. The effect of asymmetrical forces 
acting on the towed vessel, such as the presence of a 
propellor or an environmental force, are introduced through 
a bias in the data file. All dimensions are 
nondimens ionalized with respect to the towed ship's length 
between perpendiculars (LBP) . 

Three vessels were studied (Fig. 2): 

1. a 191 foot barge with a skeg, with no propellor 
(i.e. , no bias) ; 

2. a 1066 foot tanker with no skeg, but with a propellor 
(i.e., with a bias); 

3. the same barge as in 1), but without the skeg and 
with a propellor (i.e., with a bias). 

Unlike previous studies, this work includes the effect 
of athwartship position of the towline attachment point in 
the stability of the towing system. 

Chapter II provides background into the problem 
formulation and stability analysis used in this study. 
Chapter III presents the results of the analysis and 
discusses some practical aspects of these results. Chapter 
IV discusses the conclusions which can be made from the 
results of this study on the stability of the towing system 
and the use of the techniques used herein. 
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* Barge with no skeg, with bias 





* Tanker with bias 

Figure 2. Body Plans of Vessels Studied 
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II. PROBLEM FORMULATION AND METHOD OF APPROACH 



Slow motions of a towed vessel in the horizontal plane 
are described by a system of six nonlinear, coupled, 
differential equations. [Ref. 4 and 5] In its standard form 
this system is 



1 

*1 = [F 1 (x 1 ,x 2 ,x 3 ) + T s e (x 4 ,x 5 ,x 6 )] 

m " X u 
. I z - N r 

x 2 = [F 2 ( x i,X 2 , x 3 ) + T sway (x 4 ,X 5 ,X 6 ) ] 



+ -~-[F 3 (x 1 ,x 2 ,x 3 ) + x pT sway (x 4 ,x 5 ,x 6 )-y p T surge (x 4 ,x 5 ,x 6 ) ], 



Nv 



X 3 = [F 2 (x 1 ,x 2 ,x 3 ) + x p T sway (x 4 ,x 5 ,x 6 ) ] 



D 



m - 



[ F 3 (x 1 , x 2 , x 3 ) +XpT sway (x 4 / x 5 / x g) Yp^surge ( x 4 ' x 5 ' x 6^ ^ ' 



x 4 = x-^cosxg - x 2 sinx 6 - U, 



x 5 = x 1 sinx 6 + x 2 cosx 6 , 
x 6 = x 3 ' 
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where 



T surge( x 4' x 5' x 6) = R x ( x 4 ' x 5 ' x 6 ) cosx 6 +R y (x 4 , x 5 , x 6 ) sinx 6 , 
_T sway( x 4' x 5' x 6) = R x ( x 4 ' x 5 , x 6 ) sinx 6 ~R y (x 4 , x 5 , x 6 ) cosx 6 , 
D denotes the known quantity 



D = (m-Y v ) (I z -N r )-Y r N v , 



and 



F 1 ( Xl , X 2 , X 3 ) =X u X 1 +^X uu Xi+1/6X uuu Xi+^X vv X2+^X vvu X2X 1 



2 2 

+^X rr X3+^X rru x 3 x 1 + (X vr +m) x 2 X3+X rvu x 1 x 2 X3 , 



F 2 ( X 1 ' x 2 ' x 3 ) =Y 0 +Y 0u x l +Y ouu x l +Y v x 2 +1 /6Y vvv X2+^Y vrr X2X3 

+ Y vu x i x 2 + ^ y vuu x 2 x 1 + ( Y r- mx l> x 3 +1 /6Y rrr x 3 +^Y rvv x 3 x 2 
+Y ru x 3 x l + ^ Y ruu x 3 x l ' 



2 3 2 

F 3 ( X 1 ' x 2 ' x 3 ) =R 0 +R 0 u x 1 +R ouu x 1 +R v x 2 + ^/ ®vw x 2 + ^vrr x 2 x 3 

2 3 2 

+N vu x l x 2 + ^ N vuu x 2 x l +N r x 3 +1 / 6N rrr x 3 + ^ N rw x 3 x 2 

+N ru x 3 x l + ^ N ruu x 3 x l* 

In the above equations, x 1 denotes the sway velocity in 
surge (longitudinal motion) of the towed vessel, x 2 the 
velocity in sway (lateral motion) , x 3 the angular velocity 
in yaw (turning motion about the vertical axis) , x 4 and x 5 
the coordinates of the center of gravity of the towed vessel 
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with respect to an (x, y) -coordinate system moving with the 
towing vessel, and x 6 the towed vessel yaw angle. Further, U 
is the steady towing vessel velocity in the x-direction, xp 
and yp are the coordinates of the towline connection point 
on the towed vessel with respect to an (X, Y) -coordinate 
system with its origin at the towed vessel center of 
gravity, and R x , Ry are the towline restoring forces. The 
towing system configuration and notation conventions are 
shown in Figure 1. Expressions for F-^, F 2 , F 3 are derived by 
Taylor expansion in terms of the relative velocities x^, x 2 , 
x 3 of the towed vessel with respect to the water. In 
nonlinear analysis terms up to third order are used whereas 
terms beyond third order and second- and higher-order 
acceleration terms are usually neglected. Subscripts u, v, r 
indicate derivative of force-moment component with respect 
to x 3 , x 2 , x 3 respectively, and subscript c indicates 
propellor dependent terms., which represent a source of 
system asymmetry. These terms are zero in the absence of a 
propellor. Terms X a j 3c , Y a k c , N at)C , where a, b, c are dummy 
independent variables representing u, v, r, are usually 
called slow motion derivatives. In unsteady reference 
motion, slow-motion derivatives are considered as functions 
of the frequency of motion. In our study of slowly varying 
reference motions, we assume that slow motion derivatives 
are time independent. This is a good approximation for ships 
with usual hull shapes and moderate speeds. 
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R x and Ry denote restoring forces from the towline, and 
for a quasistatic towline response they are expressed as 
implicit functions of x 4 , x 5 , x g . In this study, the model 
used for the towline is that of an inextensible heavy 
caternary with nonlinear force-displacement characteristics 
as given in [Ref. 3], 

In compact notation the above system of six ordinary 
differential equations is denoted as 

x=f(x) (1) 

where x and f are six dimensional vectors. To analyze the 
stability properties of (1), the first step is to identify 
the equilibrium configuration of the system. For this we 
have to solve a system of six nonlinear, coupled algebraic 
equations 

f(x)=0 (2) 

where x denotes an equilibrium configuration. It can be 
shown [Ref. 5] that system (2) has at most three solutions 
in x corresponding to three distinct equilibrium positions. 
In this study we concentrated our efforts on one of these 
equilibrium positions, namely the one which, in the absence 
of a bias in the system, corresponds to the towed vessel 
being located directly astern of the tow- tug. This is the 
most interesting in applications. Having computed x, its 
stability properties can be established as follows: 
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Linearization of (1) around x leads to the linear system 

z=Az (3) 

where z represents the excursion from the equilibrium x, and 
A is a constant 6x6 matrix. If all eigenvalues of A have 
negative real parts, then x is stable; otherwise it is 
unstable. 

In this study, we performed parametric analysis of the 
central equilibrium in therm of towline length Lw, and the 
towing point coordinates with respect to the center of 
gravity of the towed vessel, xp and yp. These parameters can 
be easily changed before or during towing operations and can 
provide a means of passive control of the towing system. 
Parameter Lw directly affects the amount of tension 
developed by the towine. Parameter xp is directly related to 
the towline restoring force and moment. A small value for xp 
may not be able to provide adequate restoring moment and may 
not guarentee system stability. On the other hand a very 
large value for xp may result in over-compensation and 
therefore instability. Nonzero yp values result in a source 
of asymmetry introduced in the system. For a biased system 
(for example due to the presence of a propellor or 
environmental forces) , it should be expected that an extra 
appropriate bias introduced via a nonzero yp helps 
counteract the effect of the former bias, and hence improve 
stability. 
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The particular equilibrium position will lose its 
stability when an eigenvalue of the A matrix in (3) changes 
its sign from negative real part to positive real part. The 
case when a real eigenvalue crosses zero has ben analyzed in 
detail in [Ref. 5]. This corresponds to a static loss of 
stability with generation of additional equilibrium 
positions in the form of solution branching. In this study 
we focussed our attention on the case when a complex 
conjugate pair crosses the imaginary axis. This corresponds 
to a Hopf bifurcation: the particular equilibrium 
experiences a dynamic loss of stability and the system 
begins to oscillate. The resulting periodic solutions can be 
stable or unstable, but at any rate, such a situation is 
hazardous and should be avoided during towing operations. 
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III. RESULTS AND DISCUSSION 



A. BARGE WITH 8 KEG 

The first vessel to be studied was an unpowered barge 
with a skeg aft. Since there is no propellor to introduce 
a bias, the barge has athwartship symmetry. 

1. Figure 3: Critical Real Part vs. xp 

Program T0WBIF1 calculates eigenvalues for specific 
Lw and yp, creates a file for each of six real and six 
imaginary parts, and creates a separate file containing the 
largest real part. The real part with the largest value is 
the critical indicator of the system's stability: if it is 
greater than zero, the system will be unstable; if less 
than zero, the system will be stable. 

Figure 3 shows plots for the critical parts for 
yp=0.05 and three values for Lw. The region where the plot 
is greater than zero indicate that range of xp where the 
system is unstable. For example, for Lw=0.5 the critical 
real part is greater than zero for the range of xp=0.18 to 
xp=0.48, so the system is unstable within this range. 

Note that as Lw increases, the unstable range 
becomes smaller, until for Lw=3 . 0 there is no region 
greater than zero. Therefore, the system will be stable 
for all values of xp; i.e, the barge should exhibit no 
unstable motion. 
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Figure 3. Critical Real Part vs. xp - Barge with skeg 



2. Figure 4: yp vs. xp, Lw as Parameter 

Program TOWBIF2 does the same calculations as 
T0WBIF1 over a range of values of yp with a given Lw, 
instead of a single value of yp and Lw. In essence. Figure 
3 represents a cut of Figure 4 at a single value of yp and 
Lw. Unlike T0WBIF1, T0WBIF2 writes a point only where the 
critical real part changes sign. When plotted, these 
produce a curve delineating stable and unstable regions of 
the parameter space. This is the point of the process; we 
are more interested in finding what parameters produce 
stable or unstable system than the actual results of the 
equations of motion. 

Recalling Figure 3, The area inside the curves 
represent the unstable region. For example, for Lw=0.5, 
the system is unstable for all values of yp within the 
range xp=0.2 to xp=0.5. Increasing Lw first decreases the 
unstable range of xp for high yp, then decreases the 
unstable range of yp. For large Lw (>4.0), the unstable 
range virtually disappears. 

3. Figure 5: Lw vs. xp, yp as parameter 

Program T0WBIF3 performs the same calculations as 
T0WBIF2, but with Lw as the ordinate and yp as the 
parameter. Thus Figure 5 provides the same information as 
Figure 4 but with a different perspective. 
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Figure 4. yp vs. xp, Lw as parameter 
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Figure 5. Lw vs. xp, yp as parameter 



As in Figure 4, the unstable region is inside the 



curves. It clearly shows how increasing Lw decreases 
the unstable range for a constant yp, as was evident in 
Figure 4. It also shows that, for constant Lw greater than 
about 0.7, increasing yp also decreases the unstable 
region. In the narrow range of Lw from 0.2 to 0.7, 
increasing yp increases the extent of the unstable range of 
xp. This effect is apparent in Figure 4, but more 
dramatically presented in Figure 5. It would appear that 
using two views of the data would emphasize aspects of the 
curve that may be overlooked with one view. 

Since positive values of yp represent port side 
placement of the towline attachment point, and negative 
values starboard side placement, both positive and negative 
values for yp were studied. As expected from the port- 
starboard symmetry of the barge, curves for positive and 
negative values of yp were identical, and only positive 
values were presented here. 

From an operational point of view, one may conclude 
from these curves that for the unpowered barge, placing the 
towline on an attachment point to either side, as far 
forward as possible, will make the tow stable for the 
greatest range of towline length, but the towline should be 
kept no shorter than the length of the barge. 
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B . TANKER 



The second vessel studied was a tanker typical of those 
now in service. The effect of the tanker's propellor makes 
the hull asymmetrical; this effect is represented by a bias 
included in the tanker data file. 

1. Figure 6: Critical Real part vs xp 

Figure 6 plots data generated from T0WBIF1 with 
yp=0.10 and two values for Lw. These plots show the stable 
region to be between the two zero values for the curve. 

Note that the stable region becomes smaller with increasing 
Lw. 

Note also that both curves are discontinuous in their 
slopes. The critical real parts file is a composite of 
several results files, each of which is critical over a 
certain range. Each results file forms a smooth curve; 
thus the curves plotted on each T0WBIF1 figure may be 
combinations of the critical section of several results 
files. 

Finally, note that the stable region occurs over a 
narrow range of xp, unlike the barge with skeg discussed 
earlier. 

2. Figure 7: Lw vs. xp, yp as Parameter 

Figure 7 was produced from data generated by T0WBIF3 
for positive values of yp. As was shown in Figure 6, the 
stable region is inside the curves. The vertical line at 
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Figure 6. Critical Real Part vs. xp - Tanker 
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Figure 7. Lw vs. xp, yp>0 as parameter 



xp=0.404 is a common crossing point for all curves. Each 
curve is formed by two cusps, with the upper cusp 
dominating with decreasing yp. Each cusp is the plot of 
different critical pair of eigenvalues; the "nose" in the 
curves is the point where they intersect. 

Note how the stable region gets smaller with 
decreasing yp for Lw less than 0.7, for example, with 
yp=0.0 and Lw greater than 1.0, there is a very narrow 
range of xp where stability can be assured. 

3. Figure 8: Critical Real Fart vs. xp 

T0WBIF1 was again used to form Figure 8, this time 
with one value for Lw (Lw=0.6) and three negative values 
for yp. The negative yp curves pass from stable to unstable 
regions, with the stable ranges for xp getting smaller as 
xp becomes more negative. 

As in Figure 6, the curves are composites of those 
results curves which are critical over a particular range 
of xp. 

4. Figure 9; Lw vs. xp. Negative Values of yp as 
Parameters. 

Figure 9 data was generated from T0WBIF3, with yp=0.0 
curve included to provide continuity with Figure 7 . 

The stable region gets smaller as yp decreases from 
0.0. At yp=-0.10 the "nose" between upper and lower cusps 
appears to be tipping up, with the region inside the "nose" 
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Figure 8. Critical Real Part vs. xp - yp<0 
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Figure 9. Lw vs. xp, yp<0 as parameter 



being stable. The curves plotted in Figure 8 were formed 
using a value of Lw which cut through this nose, thus 
forming the sinuous curves which pass in and out of the 
stable region. Note that for the most negative values of 
yp, the cusps have disappeared, and the stable range of xp 
is slightly increasing. 

5. Figure 10: Lw vs xp. Negative Values of yp as 
Parameters. 

Figure 10 is a "close-up" of Figure 9, focusing on 
what is happening around yp=-0.10. The lower cusp tips up 
and merges with the upper to form a single curve. Note how 
the stable range of xp virtually disappears for Lw greater 
than 0.5 for yp=-0.10 and -0.11. As was seen in Figure 9, 
stable range for xp for Lw greater than 0.5 reappears with 
yp less than -0.15. 

C. BARGE WITHOUT SKEG 

The third vessel was a self-propelled version of the 
barge studied in Section A (not under power during tow) , 
but without the skeg. As with the tanker, the presence of 
the propellor, simulated by a bias in the data file, 
introduces port-starboard asymmetry. 

1. Figure ll: Critical Real Part vs xp for Lw=l.5 
Figure 11 shows curves for three values of yp 
(greater than zero, zero, and less than zero) and one value 
of Lw. The curves show the stable range of xp to be between 
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Figure 10. Lw vs. xp, yp<0, Close-up 
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Figure 11. Critical Real Part vs. xp - Barge w/o skeg 



the zero crossing points of the critical real parts, as in 
the tanker case. Also similarly to the tanker, the stable 
region increases with increasing yp. These results are 
opposite to the propellor-less barge with the skeg. 

2 . Figure 12 : Lw vs xp 

Figure 12 dramatically shows how decreasing yp 
reduces the stable region. The vertical line at xp=0.14 was 
common to all values of yp greater than and equal to zero. 
For values of yp less than zero, the smooth shape of the 
curve is apparent. 

The tanker and self-propelled barge cases 
dramatically demonstrate the effect that a bias, like a 
propellor, can introduce to the stability of the system. 

D. PRACTICAL OBSERVATIONS 

Analysis of the graphs suggests some general principles 
which may be applied when conducting slow speed towing 
operations with the vessels discussed in this chapter. 

While these principles are of course not generally 
applicable to all vessels, they illustrate how the analysis 
techniques employed in this work can be applied to other 
vessels . 

For the unpowered, symmetric barge with a skeg, the 
operator should have the towline attachment point as far 
out to either side as possible and the towline as long as 
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Figure 12. Lw vs. xp, yp as parameter 



practical. The attachment point can then be placed at any 
location forward of the center of gravity with stability 
assured. Conversely, if the attachment point must be on 
the centerline, placing it as far forward as possible 
(about half the barge's length forward of the center of 
gravity) will assure stability for all towline lengths. 

For vessels with an asymmetrical bias (e.g., with a 
propellor) , but without skegs, the attachment point needs 
to be as far to the biased side as possible (in the cases 
of the tanker and self-propelled barge, the +yp or port 
side) and placed forward of the center of gravity the 
distance indicated on the graph for all towline lengths. 
Placement of the attachment point on the opposite side (in 
the cases studied, the starboard side) will virtually 
assure the system to be unstable. 
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IV. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

This study highlighted the effect of athwartship 
position of the towline attachment point. The common 
assumption among ship operators prior to this reseach held 
that placing the towline on the centerline on the foremost 
point of the towed vessel would create the most 
advantageous towing situation. Studies such as [Ref. 2] 
have shown that towing stability can be dependent on the 
longitudinal placement of the attachment point. This 
research has shown that for certain conditions, attaching 
the towline off the centerline can also improve towed 
stability. The optimum towing configuration requires a 
combination of all three parameters - longitudinal and 
athwartship placement of the towline attachment point, and 
towline length. 

The bifurcation technique used in this study can be 
used to produce stability information useful to ship 
designers and towing operators. Stability information can 
be assembled into a convenient graphical form that clearly 
defines the regions of stable and unstable operation based 
on the parameters the operator has the most control over - 
the placement and length of the towline. 
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For the ship designer, this technique can be useful in 
determining the implications particular design decisions 
would have on the vessel's performance under tow. 

Depending on the vessel's use, adjustments to the design 
can be made to improve towing stability, or the customer 
can be forewarned to avoid certain kinds of operations. 
Since nearly all vessels are towed at some time, towing 
performance should be analyzed for all vessels. 

For the towing operator, this technique can provide 
readily available information about how a particular vessel 
will respond under tow. The operator can then adjust the 
towing parameters (i.e., placement the attachment point 
and/or length of the towline) so the tow will be in its 
most stable condition, or, if unavoidable, know that a 
particular towing situation will be potentially dangerous 
and make preparations to deal with it. 

Since ship data is inputted through a data file, the 
towed performance of any vessel can be analyzed with this 
method, including structures such as offshore oil 
platforms. Existing vessels can be analyzed, as well as 
different loading conditions. 
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Two principal disadvantages are associated with this 
technique: 

1. The programs are dependent on the quality of the 
data provided. Determining hydrodynamic coefficients 
and resistance data requires tow tank experiments and 
analysis, and are not obtained for most vessels; 

2 . The programs require large amounts of computer time 
and memory to run, which may not be available or too 
costly for potential users, especially to run 
extensive "what-if" scenarios. This problem may be 
alleviated as more inexpensive, high speed, high 
capacity micro- and personal computers become 
available. 



B. RECOMMENDATIONS 

This study was done for only one set of conditions. 
Further research can be done in determining the effect of 
varying conditions, such as different speeds or maneuvering 
by the towing vessel, on the stability of the tow. External 
forces are modelled by the bias in the data file. A 
systematized method of introducing biases into the data 
would enable the analysis of the effect of environmental 
conditions on the towing system. 

Further work should be conducted to improve the "user- 
friendliness" of the programs. As currently configured, the 
programs must be run instructively, and graphics produced 
offline. This is a time consuming process which does not 
use the full capabilities of either the programs or the 
graphics capabilities of the mainframe. Program 
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improvements should focus on streamlining computations and 
user interaction, and incorporating graphics, with the goal 
of making it available as a ship design tool. 
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APPENDIX 



Driver programs used in this thesis are shown here 
Subroutines can be obtained by contacting: 

Prof. F.A. Papoulias 

Department of Mechanical Engineering 
Naval Postgraduate School 
Monterey, CA 93943 



FILE: T0WB1FI FORTRAN A1 



C 

C 

C 

C 

C 

C 

c 

c 



c 



c 



c 



c 



PROGRAM TOWBII 

BIFURCATION ANALYSIS OF TOWING SYSTEMS 
PARAMETER DEPENDS ON IPAR 
IPAR = 1 ; XP 
2 : YP 
2 : LW 

IT NEEDS SUBROUTINES FROM TOWING. FTN 
IMPLICIT DOUELE PRECISION (A-H.O-Z) 

DOUBLE PRECISION MASSP . NVD . NV . NRD . 1 ZZ • NR . LB . LEN ,LW , 

1 NO . NOU . NOUU • NVVV . NVRR . NVDD . NVU , NVUU . NRRR . NR VV . 

Z NRDD . NRU . NR UU . ND . NDDD » NDVV . NDRR . NDU . NDUU . NVRD 

DIMENSION I V 1 ( 6 ) • A ( 6 . 6 1 . V V t 3 ) » X ( 6 ) . WR ( 6 ) . W I ( 6 ) • Z ( 6 • 6 1 »SV2 ( 6 ) 

COMMON/ 1 NTGR/ ISKEG . NREDP ♦ I TYS * ID. I FDS . I STAB . IPROP 
COMMON/SPAR /MASSP . LW » XPP . YPP • LB 
CCMM0N/SURGE/SU(7 ) 

common/xsurg/xu.xuu.xuuu 

CCMMON/SWAY/SW( 1 5 ) 

CCMMCN/YAW/ YA ( 16 ) 

COMMON/MTER/VCAR .RHO.ABS.CON1 .CON2 
COMMON/RESIST/VEL(40) ,RESI KO) 

COMMON/VELE/UEL ( 100 ) 

COMMCN/POSTN/X1 .YI .21 
COMMCN/GEOM/AL . R W . G . AET . HW . HW 1 
COMMON/PROP/ALE .P .EY.DIA.ANIU 
COMVCN/CTNR/XCC99 ) .YCC99 1-ZC(99 1 .TCl Q 9) 

COMMON/ 1 NT 1 / IC 
COMMON/DOC/UC. ALPHA 

common/uept/rlx.rly.rlz 

COMMON/SLOPE /PDRXX.PDRXY.PDRYX.PDRYY 
COMMON/S LAN /R XX 6 . RYY6 . RXX • RYY 

OPEN (UNIT=S5. FILE* ‘BARGED* .STATUS* ' OLD’ ) 

OPEN (UNIT=1 ,FILE=’RES0’ . STATUS* * NEW* ) 

OPEN ( UN I T = 1 1 .FILE*' RESIR ' .STATUS** NEW ) 

C?EN (UNIT=I2.FILE=*RESDR' .STATUS* ’ NEW' ) 

OPEN (UNIT=I3.FILE=’RES:R' ,$TATUS='NEW ) 

OPEN (UNIT=1< .FILE* ' RES4R ’ .STATUS* ' NEW ) 

OPEN C UN IT=15.FILE=‘RES5R’ .STATUS* 'NEW ) 

OPEN (UNIT=I6.FILE*‘RES6R* . STAT US* ' NEW ' ) 



TOWOOO 10 
TOWOOODO 
TOW00030 
TOWOOOAO 
TOW00050 
TOW00060 
TOW00070 
TOWOOC8 0 
TOW00090 
TOWOO 1 00 
TOWOOllO 
TOWOO 120 
TOWOO 130 
TOW00I40 
TOWOO 150 
TOWOO 160 
TOW00I70 
TOWOO 180 
TOWOO 1 9 0 
TOWOC200 
TOWOODIO 
TOW00220 
TOW00230 
TOW002AQ 
TOW00250 
TOW00260 
TOW00270 
TOW00280 
TOW00290 
TOWOOJOO 
TOWOO 2 I 0 
TOW00320 
TOW00230 
TOW00340 
TOW00350 
TOW00360 
TOW00270 
TOW00280 
TOWOO 3 90 
TOW00400 
TOW00M0 
TOWOOADO 
TOWOOA30 
TOWC0440 
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OPEN C UN I T = 2 1 . F ILE= ' RES1 1 * . STATUS* * NEW* ) 


TCW00450 

TCW00460 




OPEN (UNIT = 22 . F I LE= * RES2 I ' . STATUS 5 * NEW ) 


TOWOOATO 




OPEN (UNIT 5 23 .FILE=*RES3I * .STATUS 5 ’ NEW ) 


TOWQ0A80 




OPEN (UNIT 5 24.FILE=*RE$4I 1 .STATUS 5 * NEW ) 


TCWOOA90 




OPEN ( UN I T = 25 .FILE 5 * RES5 I ' .STATUS 5 ‘NEW ) 


TCW00500 




OPEN (UNIT=26 .FILE=*RES6I * .STATUS 5 ‘NEW ) 


TCW0C5 1 0 




CALL INPUTD(IO) 


TOWO0520 
TOWQ05 3 0 




VCAR sVCAR-l .689D0 


TCW005AQ 




MAT? =0 


TQW00550 




I FLOW= I 


TOW005 60 




WRITE (",1001) 


TOW00570 

TCW00580 




READ (",") IPAR 


TCW00590 




WRITE (".1002) 


TOW00600 




READ <".") A1.A2 


TCW00610 




WRITE (".1003) 


TOW00620 




READ (",") NUMI 


TOW0O63O 




WRITE (",1005) 


TCW00640 




READ (",") IKB 


TOW00650 




WRITE (".1006) 


TCW00660 




READ (",") NEOL 


TCWC0670 




IF ( IKB.GT.NECL) GO TO 500 


TCW00680 




DO 1 1 = 1. NUN 1 


TCW00690 




WRITE (-.2001) I, NUMI 


TCW0070Q 




AA = A 1 ♦ ( A2- A 1 )" ( 1-1 )/ (NUMI- I ) 


TOW0C710 




IF (IPAR.EO.l) XPP =AA 


TCW00720 




IF (IPAR.E0.2) YPP=AA 


TOW00730 




IF ( I PAR . EO. 3 ) LW =AA 


TOW007A0 




AL 5 LW*LB*0 . 3048 DO 


TCW0075C 




ALE 5 AL 


T CW 00760 




CALL STABIL(IVV.W.ISOL) 


TCW00770 




SET V=VV(K) FOR BIFURCATION ANALYSIS OF K-TH EQUILIBRIUM 


TCW0Q780 

TCWC0790 




IF (IW.NE.NEOL) GO TO 1 


TOW008CO 

TCW00810 




V=VV( IKB) 


TCWC0820 




IF (DABS(V).GT.I .DO) STOP 1111 


TCW00830 




CALL EOUILB (V.X.RES-RX-RY) 


TOW008AO 




CALL LINEAR (X.RES.A.RX.RY) 


TCW00850 




CALL RG(6 -6 .A.WR.WI . MAT2 . 2 . I Vl . $V2 . I ER 1 ) 


TOW00860 




IF (IER1.NE.0) STOP 2222 


TOW008 70 




CALL DEGSTB(DEOS.WR) 


TOWO0880 




WRITE (1.10) AA.DEDS 


TOW0089C 




DO 11 J=I .6 


TCW00900 




JR= 1 0 ♦ J 


TOW00910 




WRITE (JR. 10 ) AA.WR(J) 


TCW00920 




JI 5 20*J 


TOW00930 




WRITE (JI.10) AA.WI(J) 


TOW009A0 


1 1 


CONTINUE 


TOWO0950 


1 


CONTINUE 


TCWO0960 


500 


STOP 


TOW00970 


10 


FORMAT (2D20.IO) 


TCW00980 


1001 


FORMAT C ENTER 1 : XP VARIATION 1 ,/, 


TCW00990 




1 * 2 : YP VARIATION 1 ./, 


TCWO 1000 




2 1 3 : LW VARIATION 1 ) 


TOWQIOIO 


1002 


FORMAT C ENTER PARAMETER RANGE 1 ) 


TOW01020 


1003 


FORMAT C ENTER NUMBER OF INCREMENTS 1 ) 


TOWO 1030 


1005 


FORMAT (‘ ENTER EOUILIERIUM NUMBER*) 


TOWO 1040 


1006 


FORMAT C ENTER ESTIMATED NO. OF EQUILIBRIA*) 


TOWO 1 050 


2001 


FORMAT (215) 


TOWO 1060 




END 


TOWO 1070 
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T0WB1F2 FORTRAN A1 



FILE: 



C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c 



c 



c 



c 



c 



PROGRAM TOWBI2 
PROGRAM TOWBIF.FTN 

BIFURCATION ANALYSIS OF TOWING SYSTEMS 

PARAMETERS ARE: Xp. Yp 

IT NE-EDS SUBROUTINES FROM TOWING. FTN 

USER DEPENDENT SUBROUTINES: 

DEGSTB = CURVES ENCLOSING REGION II OF FIGURE 13 
(SUBROUTINE DEGSTB IS IN TOWING. FTN) 

DS1 ^ CURVES ENCLOSING REGION V OF FIGURE 13 
(SUBROUTINE DS1 IS IN SPMB1F) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DOUBLE PRECISION MASSP . NVD . NV , NRD . 1 ZZ , NR , LB .LEN ,LW . 

1 NO • NCU • NOUU . NVVV . NVRR . NVDD • NVU . NVUU , NRRR , NR VV . 

2 NRDD . NRU . NRUU » ND . NDDD • NDVV • NDRR .NDU.NDUU.NVRD 

DIMENSION IV1(6).A(6.6).VV(3).X(6).WR(6).WI(6)»Z(6.6). SV2 ( 6 ) 

COMMON/ 1 NTGR/ I SKEG . NREDP > ITYS > ID » IFD3 > I STAB . I PROP 
COMMON/SPAR /MASSP . LW . XPP ,YPP .LB 
CCMMON/SURGE/SU(7) 

COMKCN/XSURG/XU »XUU» XUUU 
CCMMCN/SWAY/SW( 15 ) 

COMMON/YAW/YA ( 16 ) 

COMMON/MTER/VCAR • RHO . AES . CON 1 ,CON2 
COMMON/RES IST/VEL (AO) .RES I (AO) 

COMMON/VELE/UELUOO) 

CCMMCN/POSTN/XI , Y1 ,Z1 
COMMCN/GEOM/AL • RW . G • AET .HW.HW1 
COMMON /PROP/ALE .P .EY.D1A.ANIU 
COMMON ^CTNR/XC ( 99) • YC ( 99 ) • ZC ( 9 9 ) * TC ( 9 9 ) 

COMMON/ I NT 1 / I C 
COMMGN/DOC/UC. ALPHA 
COMMON/ UEPT /RLX . RLY .RLZ 

commcn/slope/pdrxx.pdrxy.pdryx.pdryy 

COMMCN/SLAN/RXX6 • RYY6 • RXX . R YY 



TOWOOOIO 
TOWOOOZC 
TOW0003C 
TOWOOOAO 
TOWOOOBO 
TOW00060 
TOW00070 
TOWOOOBO 
TOW00090 
TOWOOIOO 
TOWOOllO 
TOW00120 
TOW00130 
TOWOOIAO 
TOW00150 
TOW00160 
TOW00170 
TOW00180 
TOW00190 
TOWOOZOO 
TOWOOZIO 
TOWOOZZO 
TOW00220 
TCW002A(j 
TOW00250 
TOW00260 
TOW002 7 0 
TOW00280 
TOW0C2 90 
TOW00300 
TOW00210 
TOW00220 
TOW002 20 
TOW003A0 
TOW00250 
TOWOO 260 
TOW00270 
TOWOO 28 0 
TQW00290 



OPEN 


( UN 1 T = 25 


• FILE** BGKEG2 * »STATU3= ’OLD' ) 


TOWOOAOO 


OPEN 


( UNI T= I 1 


• F 1LE= ’ PE31R * »STATUS= ' NEW* ) 


TOWOOAIO 


OPEN 


( UN I T = 1 2 


• FILE= 'RE32R 1 » STATUS= * NEW ' ) 


TCW00A20 


OPEN 


( UN I T = 1 3 


• FlLEs’RESSR’ • ST ATUS = * NEW ’ ) 


TOW00A20 


OPEN 


(UNIT = 1A .FILE=’RESAR‘ * ST ATUS= ’ NEW ’ ) 


TOWOOAAO 


CPEN 


( UN1T = 15 


• P 1 LE= ‘ RE35R ’ . STATUS* ' NEW 1 ) 


TOWOOA50 


OPEN 


( UNl T= 1 6 


• FILE* ’RES6R * , STATUS =' NEW ’ ) 


TOW00A60 








TCW00A70 


CALL 


1 NPUTD (10) 


TOWOOA0O 


VCAR 


= VCAR * 1 . 


689D0 


TOW00A90 


AL 


sLW-LB-0 


. 30A8D0 


TOW00500 


ale 


= AL 




TOW005I0 


MATZ 


= 0 




TOW00520 


I FLOW 


= 1 




TOW005 2 0 


EPS 


=1 .D-5 




TOW005A0 


1LMAX 


= 1500 




TOW00550 








TOW00560. 


WRITE 


<•*, 1001 ) 


TCW00570 


READ 


("••<) 


A 1 , A2 


TOW0058 0 


WRITE 


(•*,1002 


) 


TOW00590 


READ 


(«,*•) 


NUM1 


TOW00600 


WRITE 


(•*.1002 


) 


TOWOO k > 1 0 


READ 


(•*.*•) 


B 1 . B2 


TOW0C620 


WRITE 


(« . 100A ) 


TCW00630 


READ 


(**,**) 


num: 


TOWOOeAO 
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WRITE (-.1005) 

READ (".*) 1KB 

WRITE C. I 00 61 
READ (",") I DC 

DO 1 1=1 ,NUM1 

WRITE (**,2001 ) l.NUMI 

VP P = A 1 * ( AC- A I )M1-1 l/CNUMl-l ) 

XPP =B 1 

CALL STABILC1VV.VV.1S0L) 



V = VV( 1KB) 

IF CDAEGCV) .GT.1.D0) STOP 1111 
CALL EOUIL8(V.X.RES.RX,RY) 

CALL LINEAR (X.RE3.A.RX.RY) 

CALL RGC6 .6 .A.WR.Wl .MATZ -Z. 1V1 .SV2 

IF (IER1.NE.0) STOP Z2Z2 

IF (1DS.E0.1) CALL DEGSTB(DEOS,WR ) 

IF (IDS. EG. 2) CALL D31CDE0S.WR) 

DE0GNN=DE03 

XPNN=XPP 

PR=DEOCOO*DEOGNN 
IF (PR.GT.O.DO) GO 70 3 
L = L*1 

IF (L.GT.6) STOP 1000 
1L = 0 

XPO=XPOO 

XPN=XPNN 

D£OGO=DEOSOO 

DEOGN=DEOGNN 

XPL=XPO 

XPR=XPN 

DEOSL = DECGO 

DEOSR=DEOSN 

XPP = ( XPL * XPR )/2 . DO 

CALL STAB1LC 1VV.VV. 1S0L) 

V = VVC I KB ) 

IF CDABS(V) .GT. 1 .DO) STOP 1111 
CALL ECUILE(V.X.RES.RX.RY) 

CALL LINEAR(X.RES.A.RX.RV) 

CALL RG(6 .6.A.WR.W1 .MATZ.Z. I VI .SV2 
IF C1EP1.NE.0) STOP 2222 
CALL DEGGTB(CEOS.WR) 

DE0SM=DE03 



TOW00650 
TOW00660 
TOW00670 
TOW00A80 
TOW00690 
TOW00700 
TOW 0 0 7 1 0 
TOW00720 

TOW00730 
TOW007A0 
TOW00750 
TOWOO 7 60 
TOW00770 
TOW007B0 
TOWOO 7 90 
TOWOOBOO 
TOW00810 
TOW00820 
TOW00830 
TOWOOBAO 
TOW00850 
TOW00860 
TOW00870 
TOW00880 
TCWC0B90 
TOWOO 9 0 0 
TOWOO 9 1 0 
TCW00920 
TOW009 30 
TOW009A0 
TOW00950 
TOW00960 
TOWOO 9 70 
TOWOO 98 0 
TOW00990 
TOWOIOOO 
TOWOIOIO 
TOWOIOZO 
TCW01030 
TCWOIOAO 
TOW01050 
TOW01060 
TOW01070 
TOWOIOBO 
TOW01090 
TOW01 100 
TOWOl 1 10 
TOW01120 
TOWOl 130 
TOWOllAO 
TOWOl ISO 
TOWOl 160 
TOW01170 
TOW01180 
TOW0I190 
TCW01200 
TOW01210 
70W01220 
TOW01230 
TOW012A0 

* I ER 1 ) TOW01250 

TOW01260 
TOWOl 270 
TOW01280 



SET V = VV ( K ) FOR BIFURCATION ANALYSIS OF K-TH EQUILIBRIUM 
V = VV( IK8) 

IF (DABS(V).GT.l.DO) STOP 1111 
CALL EOUILB(V.X.RES.RX.RY) 

CALL LINEAR(X.RES.A.RX.RY) 

CALL RGC6.6.A.WR.W1 .MATZ.Z. I VI .SV2.IER1 ) 

IF (IER1.NE.0) STOP 2222 
IF (1DS.EQ.1) CALL DEGST8 ( DECS . WR ) 

IF (IDS.E0.2) CALL DSKDEOS.WR) 

DECGCO=DECS 
XPOO =XPP 
L =0 

DO 2 J=2 . NUM2 
WRITE {".«) J 

XPP=B1*(B2-B1)"(J-1 )/(NUM2- 1 ) 

CALL 3TAEIL ( 1VV » VV. I SOL ) 

SET V = VV ( K ) FOR BIFURCATION ANALYSIS OF K-TH EQUILIBRIUM 



1ER 1 ) 



XPM=XPP 


TOW01290 


PRL-DE03L-DE0CM 


TOW01300 


PRR=DEOSR*DEOSM 


TOW01 310 


IF (PRL.GT.O.DO) GO TO 5 


TOW01320 


XPO=XPL 


TOW01320 


XPN=XPM 


TCW013A0 


DE030=DE0SL 


TOW01350 


DE0SN=DE03M 


TOW01360 


IL= IL 4 1 


TOW01370 


IF (IL.GT. ILMAX) STOP 5100 


TOW01380 


DIF=DABS(XPL-XPM) 


T OWO 1390 


IF CD1F.GT.EPS) GO TO 6 


TOW01A00 


XP=XPM 


TOW01A10 


GO TO A 


TOW01A20 


5 IF ( PRR . GT . 0 . DO ) STOP 5200 


TOW01A30 


XPO=XPM 


TOW01AA0 


XPN=XPR 


TOW01A50 


DEOSO = CEOSM 


T OWO 1 A 6 0 


DEOSN=DEOSR 


TOW01A70 


1 L = 1 L ♦ I 


TCW01A80 


IF (IL.GT. ILMAX) STOP 5100 


TOW01A90 


DIFsDABS(XPM-XPR) 


TOW01500 


IF (DIF.GT.EPS) GO TO 6 


TOW01510 


XP=XPM 


TOW01520 


A LLL = 1 0 *L 


TOW’015 30 


WRITE (LLL. 10) XP.YPP 


TOW015A0 


5 XPOO=XPNN 


TOW01550 


DEOSOO=DEOSNN 


TOW01560 


2 CONTINUE 


TOW01570 


1 CONTINUE 


TOW01580 


STOP 


TOK'O 15*0 


10 FORMAT (2D20. 10) 


TOW01600 


1001 FORMAT C‘ ENTER RANGE OF Yp VARIATION’) 


TOW01610 


1002 FORMAT (' ENTER NUMEER OF INCREMENTS IN Yp ' ) 


TOW01620 


1005 FORMAT (’ ENTER RANGE OF Xp VARIATION’) 


TOW01630 


1 00A FORMAT C ENTER NUMBER OF INCREMENTS IN Xp’) 


TCW016A0 


1005 FORMAT (’ ENTER EQUILIBRIUM NUMEER 1 ) 


TOW01650 


1006 FORMAT ( * ENTER DEGREE OF STABILITY CONTROL*) 


T OWO 16 6 0 


2001 FORMAT (215) 


TOW01670 


END 


TOW01680 


C 


TOW01690 


SUBROUTINE DSl(DEOS.WR) 


TCW01700 


IMPLICIT DOUBLE PRECISION (A-H.O-Z) 


TOW01710 


DIMENSION WR ( 6 ) 


TOW01720 


DEOS1 =- 1 . D50 


TOW01720 


DO I 1=1.6 


TOW017A0 


IF (WR( 1 ) . LT . DEOS 1 ) GO TO l 


TOW01750 


DEOGl=WR( I ) 


TOW01760 


I J = 1 


TOW01770 


I CONTINUE 


TOW0I780 


DEOC2 = - 1 . D50 


TOW01790 


DO 2 1=1.6 


TOW01800 


IF (1J.EO.I) GO TO 2 


TOW01810 


IF (WR(I ).LT.DEOS2) GO TO 2 


TOW01820 


DEOS2 = WR ( 1 ) 


TOW01830 


1 J J= I 


TOW018A0 


2 CONTINUE 


TOW01850 


DEOS=- 1 . D50 


TOW01860 


DO 2 1=1.6 


TOW01870 


IF (1 .EC. IJ.OR.I.EO.IJJ) GO TO 5 


TOW01880 


IF (WR(l).GE.DEOS) DEOS = WR ( I ) 


TOW01890 


5 CONTINUE 


TOW01900 


RETURN 


TCW0 1910 


END 


TCW01920 
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FILE 



T0WBIF3 FORTRAN A1 



C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 



c 



c 



c 



c 



c 



c 



PROGRAM TOW8IF3.FTN 

BIFURCATION ANALYSIS OF TOWING SYSTEMS 

PARAMETERS ARE: Xp . Lw 

IT NEEDS SUBROUTINES FROM TOWING. FTN 

USER DEPENDENT SUBROUTINES: 

DEGSTB = CURVES ENCLOSING REGION II OF FIGURE 13 
(SUBROUTINE DEGSTB IS IN TOWING. FTN) 

D31 = CURVES ENCLOSING REGION V OF FIGURE 13 

(SUBROUTINE DS1 IS IN SPMBIF) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DOUBLE PRECISION MASSP , NVD . NV , NRD. IZZ , NR . LB . LEN , LW , 

1 NO . NOU • NOUU • NVVV * NVRR . NVDD » NVU • NVUU • NRRR ♦ NR VV , 

2 NRDD.NRU.NRUU.ND.NDDD.NDVV.NDRR.NDU.NDUU.NVRD 

DIMENSION 1 V 1 ( 6 ) . A ( 6 . 6 ) . VV ( 3 ) . X ( 6 ) . WR ( 6 ) . W 1 ( 6 ) , Z ( 6 . 6 ) , SV2 ( 6 ) 

COMMON/ 1 NT GR/1SKEG.NREDP. ITYS. ID. IFDS. ISTAB. IPROP 
ccmmon/spar/massp .LW.XPP.YPP ,lb 
COMMON/SURGE/SUC? ) 

COMMON/XSURG/XU.XUU.XUUU 
COMMON/SWA Y/SW( 15 ) 

COMMON/ YAW/YA( 16 ) 

COMMON/MTER/VCAR .RHO.ABS.CON1 ,CON2 
CCMMON/RES IST/VEL (CO).RESI(AO) 

COMMON/VELE/UEL (100 ) 

COMMON/POSTN/XI .Y1 .Zi 
CCMMON/GEOM/AL .RW,G,AET,HW,HW1 
COMMON/ PR CP /ALE . P .EY.DIA.AN1U 
COMMON/ CT NR /XC( 99) *YC( 99 ) . ZC( 99 ) »TC( 99 ) 

CO M MCN/ 1NT 1 / 1 C 
COMMON/DOC/ UC .ALPHA 
COMMON/ UEPT/RLX . RLY . RLZ 

commcn/slope/pdrxx.pdrxy.pdryx.pdryy 

COMMON/SLAN/RXX6 .RYY6 .RXX.RYY 

OPEN ( UN1T= 35 .FILE* ’ TANKER 2 ’ . ST ATUS= ’ OLD ’ ) 

OPEN (UN1T = 1 ,F1LE= ’RE31R' *$TATUS= * NEW* ) 

OPEN (UNIT=2 .FILEs’RESDR ' .STATUS** NEW’ ) 

OPEN (UN1T = 3 .FILE=’RES3R’ . CT ATUS= ' NEW * ) 

OPEN (UN1T = A .FILEs’RESAR* , STATU3= ’ NEW * ) 

CALL 1NPUTDC10) 

VCAR = VCARM . 689D0 

AL sLW«LB"0 .30A8D0 

ALE =AL 

MATZ =0 

I F LOW s 1 

ILMAX= 1500 

EPS * 1 . D- 5 



WRITE (".1001) 

READ (".") A1.A2 

WRITE (-,1002) 

READ (",") NUM 1 

WRITE (".1003) 

READ (".") B1.E2 

WRITE (".IOOO 
READ {»,") NUM2 

WRITE (".1005) 

READ (".") IKB 



TOW0001 0 
TCW00020 
TCW00030 
TOW000A0 
TCW00050 
TOW00060 
TOW00070 
TOW00080 
TOW00090 
TOW00100 
TOWCOllO 
TOWOOIZO 
TOW00130 
TOWOOIAO 
TOW00150 
TOW00160 
TCW00170 
TCW00180 
TOW00190 
TOW00200 
TOWOOZIO 
TOW00220 
TOW00230 
TOWOOZAO 
TCW00250 
TCW00260 
TOW00270 
TCW00280 
TCW00290 
TOW00300 
TCW00310 
TOW00320 
TOW00330 
TOW003C0 
TOW00350 
TOW00360 
TCW00370 
TOWC0380 
TCW00390 
TOWOOAOO 
TOWOOA 10 
TOWOOAZO 
TOWOOASO 
TOWOOAAO 
TOWOOA50 
TOW00460 
TOWOOA70 
TOW00R80 
TOW00C90 
TOW005CO 
TOW0C510 
TOWOOSCO 
TOW00530 
TOWOOSCO 
TOW0C550 
TOW00560 
TOW0C570 
TOW005S0 
TOW0C590 
TOWC0600 
TOW00610 
TOW3C620 
TOW00630 
TCW006C0 



39 



WRITE (-1006) 



TOW0G650 

READ (-■-) IDS TOW00660 

DO I 1 s 1 » NUM I TOW006 70 

WRITE (-.2001) IiNl'Ml TOW00680 

LW = A 1 * ( A2- Al)-(I-l)/ (NUM 1 - I ) TOW00690 

AL =LW-LB-0 . 3048D0 TOW00700 

ALE=AL TOW00710 

XPf> = Bl TOWOO 720 



C 

c 

c 



V = VV( IKB ) 

IF (DABS(V).GT. 1.D0) STOP 1111 
CALL E0UIL8(V.X.RES.RX,RY) 

CALL LINEAR(X.RES.A.RX.RY) 

CALL RG(6,6.A.WR,WI . MATS • 2 • I VI . SV2 

IF ( 1EP 1 . NE . 0 ) STOP 2222 

IF ( IDS . EC. 1 ) CALL DEGST8 ( 2EOS • WR ) 

IF (IDS. EC. 2) CALL DSl(DEOS.WR) 

DEOSNN=DECS 

XPNN=XPP 

PR = DEOGOO- DEOSNN 
IF (PR.GT.O.DO) GO TO 3 
L = L* 1 

IF (L.GT.4) STOP 1000 
IL = 0 

XPO=XPOO 

XPN=XPNN 

CEOSO=DEOSOO 

DECSN=DEOSNN 

XPL = XPO 

XPR=XPN 

DEOSL=DEOSO 

DEOSR=DEOSN 

XPP= (XPL^XPR )/2 .DO 

CALL STAB I L ( I VV . VV • I SOL ) 

V = VV( IKB ) 

IF (DABS(V).GT.l .DO ) STOP 1111 
CALL E0U1L8(V-X»RES-RX.RY) 

CALL LIN£AR(X.RES.A,RX,RY) 

CALL RG(6.6.A,WR.W1 . MATS . Z . 1 VI . SV2 
IF ( 1ER 1 . NE . 0 ) STOP 2222 
CALL DEGST8(DEOS.WR ) 

DECSM=DECS 

XPM=XPP 



TOW00730 
TOW007A0 
TOW00750 
TOWOO 7 6 0 
TOWOO 7 7 0 
TOW00780 
TOW00790 
TOW00800 
TOW00810 
TOW00820 
TOW00830 
TOW00840 
TOW00850 
TOW00860 
TOW00870 
TOW00880 
TOW0C890 
TOW00900 
TOW00910 
TOW00920 
TOWOO 9 3 0 
TOW00940 
TOW00950 
TCW00960 
TOWOO 9 70 
TOW00980 
TOW00990 
TOW01000 
TOW01010 
TCW01020 
TOW01030 
TOW010A0 
TOWO 1050 
TOW01060 
TCWO1O70 
TOW01030 
TOWO 1090 
TCW01 100 
TOWO 1110 
TOW01 120 
TOW01 130 
TOWO 114 0 
TOWOl 150 
TOWO 1 160 
TOWOl 170 
TOWOl 180 
TOWO 1190 
TOWO 1200 
TOWO 1 2 1 0 
TOWO 1220 

1ER 1 ) T OWC 12 3 0 

TOWO 1 24 0 
TOWO 125 C 
TOWC126C 
TOWO 1270 



CALL STA8IL( IVV »VV . 1SOL ) 

SET V=VV(K) FOR BIFURCATION ANALYSIS OF K-TH E0UILI8R1UM 
V=VV( 1KB) 

IF (DABS(V).GT.l.DO) STOP 1111 
CALL E0UIL8(V.X.RES.RX,RY) 

CALL LINEAR (X. RES. A. RX-RY) 

CALL RG(6-6.A.WR.WI .MATZ.2. 1V1 .SV2.1ER1 ) 

IF ( 1ER 1 . NE . 0 ) STOP 2222 
IF (IDS.EO-1) CALL DEGST8 ( DECS . WR ) 

IF ( IDS. EC. 2 ) CALL DSl(DEOS.WR) 

DEOSOO=DEOS 
XPCO =XPP 
L =0 
DO 2 J = 2 . NUM2 

XPP = E1 ♦(82-81 )-(J-l )/ ( NUM2- 1 ) 

CALL STAS I L ( 1 VV » VV . I SOL ) 

SET V=VV ( K ) FOR 8 IFURCAT ICN ANALYSIS OF K-TH EQUILIBRIUM 



IER1 ) 
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PRL=DE0SL"DE0SM 


TOW01280 




PRR=DE0SR*DE0SM 


TOW01290 




IF (PRL.GT.0.D0) GO TO 5 


TOW01300 




XPO=XPL 


TCW01210 




XPN=XPM 


TCW0 1320 




DEOSO=DEOSL 


TOW01330 




DEOSN=DEOSM 


TOW01340 




1L= IL* 1 


TOW01350 




IF (IL.GT. ILMAX) STOP 3100 


TOW01360 




DIF=DABS(XPL-XPM) 


TOW01370 




IF (DIF.GT.EPS) GO TO 6 


TOW01380 




XP=XPM 


TOW01390 




GO TO 4 


TOW01400 


5 


IF (PRR.GT.0.D0) STOP 3200 


TOW01410 




XPO=XPM 


TOW01420 




XPN=XPR 


TOW0I430 




DEOGO=DEOSM 


TOW01440 




DEOSN=DEOSR 


TOW01450 




IL- IL* 1 


TOW01460 




IF ( IL.GT. ILMAX ) STOP 3100 


TOW0I470 




DIF=DABS(XPM-XPB ) 


TCW01480 




IF (DIF. GT. EPS) GO TO 6 


TOWO 1490 




XP=XPM 


TCWC1500 


4 


WRITE (L.10) XP.LW 


TOW01510 


3 


XPOO=XPNN 


TOWO 1520 




DEOSOO=DEO$NN 


TOWO 1 530 


2 


CONTINUE 


TOWO 1540 


1 


CONTINUE 


TOWO 1550 




STOP 


TOWO 1 5 60 


10 


FORMAT (2D20.10) 


TOW01570 


100 1 


FORMAT C ENTER RANGE OF Lw VARIATION*) 


TCW01580 


1002 


FOPMAT (’ ENTER NUMEER OF INCREMENTS IN Lw * ) 


TOWO 1590 


1003 


FORMAT (* ENTER RANGE OF XP VARIATION') 


TOWO 1600 


1004 


FORMAT (’ ENTER NUMBER OF INCREMENTS IN Xp') 


TCW01 61C 


lO 

o 

o 


FORMAT (’ ENTER EQUILIBRIUM NUMBER’) 


TOWO 1620 


1006 


FORMAT (' ENTER DEGREE OF STABILITY CONTROL') 


TOWO 1 6 3 0 


2 CO 1 


FORMAT (215) 


TOWO 1 640 




END 


TOWO 1 65 0 
TOWO 1660 




SUBROUTINE DS 1 ( DEOS • WR ) 


TOWO 16 70 




IMPLICIT DOUBLE PRECISION (A-H.O-Z) 


TOWO 1680 




DIMENSION WR ( 6 ) 


TCW0I690 




DEOS 1 =- 1 . DO 0 
DO 1 1=1.6 


T OWC 17 00 
TOWO 1710 




IF CWRU ) ,LT . DEOS1 ) GO TO 1 


TCW01720 




CEOS 1 =WR ( 1 ) 


TOW01730 




IJ = I 


TCW01740 


1 


CONTINUE 


TOWO 1 750 




DEOS2=-I .D30 


TCW01760 




DO 2 1=1.6 


TOWO 1770 




IF ( I J.EO. I ) GO TO 2 


TOWO 1780 




IF (WR( 1 ) .LT.DECG2 ) GO TO 2 


TOWO 1790 




DE0S2 =WR ( 1 ) 
1JJ=1 


TOWO 1800 
TCW01810 


2 


CONTINUE 


TOWO 1820 




DEOS = - 1 . D3 0 


TOWO 1830 




DO 3 1=1.6 


TOWO 1840 




IF (I.E0.1J.0R.1.E3.UJ) GO TO 3 


TOWO 1850 




IF ( WR ( 1 ) . GE . DECS ) DEOS=WR(l) 


TOWO 18 60 
TOWO 1870 


3 


CONTINUE 

RETURN 

END 


TCW0I880 
TOWO 18 9 0 
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